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NUMERICAL IMPLEMENTATION OF THE MULTISCALE AND 
AVERAGING METHODS FOR QUASI PERIODIC SYSTEMS 

TAL KACHMAN, SHMUEL FISHMAN, AND AVY SOFFER 


Abstract. We consider the problem of numerically solving the Schrodinger 
equation with a potential that is quasi periodic in space and time. We in¬ 
troduce a numerical scheme based on a newly developed multi-time scale and 
averaging technique. We demonstrate that with this novel method we can 
solve efficiently and with rigorous control of the error such an equation for 
long times. A comparison with the standard split-step method shows sub¬ 
stantial improvement in computation times, besides the controlled errors. We 
apply this method for a free particle driven by quasi-periodic potential with 
many frequencies. The new method makes it possible to evolve the Schrodinger 
equation for times much longer than was possible so far and to conclude that 
there are regimes where the energy growth stops in-spite of the driving. 


Numerical solution Time dependent potentials Multiscale averaging 

1. Introduction 

A method for the study of the dynamics for the Schrodinger equation with time 
dependent potentials p] is implemented numerically. The potential is quasiperiodic 
in both space and time. The power of this new method is demonstrated. Explo¬ 
ration of the dynamics for such potentials is motivated by experiments in optics [5] 
where hypertransport, namely transport faster then ballistic was found experimen¬ 
tally and numerically in some regimes. In theoretical work that followed 0!3i2[S] 
a classical theory was developed for the potentials relevant for these optics exper¬ 
iments. In particular, it was found in the framework of classical mechanics, that 
for smooth potentials the spreading in momentum as function of time stops. The 
calculations of the present paper are for such potentials. For short times, relevant 
for the existing experiments it was found that wave or quantum and classical dy¬ 
namics agree in general features. 

For long times it turned out impossible to compute numerically the quantum dy¬ 
namics using the standard methods 0 0, while with the method introduced in 
PJ and implemented here, calculations for such long times are feasible as will be 
demonstrated in this paper. Such calculations and comparison with the classical 
results, is of fundamental importance for the issue of quantum classical correspon¬ 
dence. The main objective of the present paper is to demonstrate the power of the 
method introduced in [T] for a physically relevant example. 

Potentials which are quasiperodic both in space and time, can manifest a high 
degree of complexity and are subject of many studies over the last decade, mostly 
in the framework of classical physics 0 0 QS HU CE2 H3 HU- With different 
experimental realizations of such potentials there is also a need for a numerical 
approach to investigate them and their asymptotic behavior. The standard way 
to numerically solve problems with time dependent potentials is based on either 


l 
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spectral methods [15j |T6] or explicit/implicit finite difference schemes [7l ITT]. In 
this paper we implement numerically a recently developed rigorously controlled 
multi-time scale averaging technique [T]. The above mentioned method has two 
distinct advantages. The first is that at each step of the averaging hierarchy there 
is a well defined and completely known bound on the numerical error. The second 
advantage and one which has far greater impact is the reduction in computational 
time accompanied with each of the hierarchical averaging steps. This reduction 
enables us to go to long time scales, impossible by the methods we are aware of. In 
this paper we describe and present an implementation for a specific problem. 

The outline of the paper is as follows. In Section 2, we present the model 
for which the method of multiple scale and averaging (MSA) is implemented. This 
model is of physical interest and importance. The MSA method that was developed 
in Ref. PQ and details of its implementation are presented in Secs. 3 and 4. The 
MSA method involves a hierarchy of computations where the first level is described 
is Sec. 3 and the following levels are presented in Section 4. The needed level is 
determined by the required precision and the time over which the system is evolved. 
In Section 5 we demonstrate that the MSA method is superior compared to the 
standard split step method, since for the same precision it is much faster. Moreover 
for the split step method there is only an empirical estimate on the error while for 
the MSA there is a rigorous bound. Using it in Section 6, we show that with the 
help of this method we could solve the Schrodinger equation for the model problem 
presented in Section 2, for an extremely long time; we conclude that for this model 
the energy does not grow to infinity in-spite of the time dependent driving potential. 
The spreading in the quantum case is wider than in the corresponding classical 
system. This results from the fact that initially we observe a lot of spreading in 
the quantum case, while the classical case does not show spreading. 


2. The model that will be studied 


In this section we introduce the model for which our MSA method developed 
in [X is implemented. In the first subsection its relation to physical systems is 
explained, while in the second one it is reduced to a form for which the MSA 
method can be applied(see Eqs. |2TlpT3t . 


2.1. The physical model. The random potentials which are prepared in optics 
mua and atom optics Hanoi mi experiments, are described by a sum of random 
Fourier components. In experiments, potentials which are composed out of a large 
number of random independent Fourier components N, are created . 

More specifically here, the Schrodinger equation for a potential 

1 N 

(2.1) V (x', t) = —= ^ A n exp ( i(k n x' - J n r )) + C.C 

V ^ n = 1 

is used. The A rn are independent, identically distributed complex random variables, 
where C.C stands for complex conjugate. The expectation values of these variables 
satisfy 

(2.2) < A m >=< A m A n >= 0 < A m A n * >= a 2 S nm 

We will study the specific model where A n = A e z ^ m with ifr m uniformly distributed 
in the interval [—7r, 7r] and A > 0. The distribution of uj n and k n is specified in 



NUMERICAL IMPLEMENTATION OF THE MULTISCALE AND AVERAGING METHODS FOR QUASI PERIODIC SYSTEMS3 


Section 5. 

The equation of motion is the time dependent Schrodinger equation 


(2.3) 
Where 

(2.4) 


id T tl> {x 1 , t) = Hvp (a/, r). 


H = ip' 2 + V{x',t), 


p' = -iS7 x >. 

In the following section we will introduce the reduction of the problem to a form 
where the MSA technique is applicable. 


2.2. Reduction of the problem. In the model with potential given by (2.1), the 


particle is expected to be accelerated most effectively in the regime of velocities 
where the Chirikov resonances 


(2.5) 


Jr) _ _n 


are formed j22 ] [23]. 

In this paper we choose Vn ' and k n uniformly distributed in the intervals , u( max l] 

and [fc( mln ), /~( max )J respectively. The crucial point is that the Chirikov resonant 
velocities are bounded in a phase space strip. This is typically the case for smooth 
potentials, and will be assumed in this paper. 

Here we would like to study what is the acceleration of a particle prepared with 
momentum or velocity v (we assume unit mass), so that all Vn ^ are far from v. 
For the classical corresponding system we found that the acceleration is negligible 
mum- Here we study the corresponding quantum mechanical system. For this 
purpose it is convenient to work in a frame of reference where the initial velocity 
of the particle vanishes. For this we perform the Galilean transformation 

(2.6) x = x' — vt, 

where v is the velocity of the moving frame (in later stages we will relate this 


quantity to the required small parameter) on the potential (2.1) 


(2.7) 


V(x- 


A N 

■ VT , T) = —= y^ COS {k n {X + VT) - J n T + (\> n ) 
V N , 


n= 1 
N 


//v 

Now, re-scaling time as 

( 2 . 8 ) 

the potential takes the form 


= — 7 = COS {knX + {k n v -u' n )T + </>„) . 


n —1 


(2.9) 

and defining 

( 2 . 10 ) 


A 


N 


vix ’ t) = 7NE 


t = TV, 


, {U)' - k n V) , 

cos | k n x - -t + i 


{u' n - k n v) 
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The potential in the moving frame takes the form 

. N 


( 2 . 11 ) 


V (X, t ) = — = C0S (k n X - Vnt + <j>n) 

V ly „_i 


The time dependent Schrodinger equation is 

(2.12) idtip (x,t) = -H(p,x,t)i[> (x,t) 

v 

Introducing the small parameter 8 = ^: 

(2.13) idt'4> (x, t) = 8H (p, x , t)ip (x, t ) 

From this point on we will use the small parameter /3 to perform the averaging 
steps introduced in the next section. In what follows we will also relate the small 
parameter p to the time scale on which we average. The Hamiltonian will be 
approximated by a finite matrix (in space and momentum) and it will be verified 
that the spreading never reaches the boundaries set by this basis. 

3. The averaging scheme 

The multiscale averaging method is based on replacing the original Hamiltonian 
by a hierarchical set of averaged Hamiltonians. In each step we perform a “peel-off” 
transformation and average a part of the Hamiltonian, for a chosen time interval 
of length T 0 = = y/v- This choice is not unique, but as shown in [T] leads to 

effective error bounds. We use the fact that Eq (2.13) is of the form of (2.1) in [T]. 


In this section the implementation of the MSA method of jTJ for equation (2.13) 
will be presented. In Appendix A we summarize the main results of Ref. [I] 

3.1. Zero order. The zero order average on the jth time interval has the form 

{V (t) = V (x,t)) 

_. l rU+i)T 6 

(3.1) V 0 ^ = - V(t)dt 


0 JjT 0 


In the case of the potential (2.1), the zero order averaging can be performed ana¬ 
lytically 

r(j+l)T 0 


Vn U) = ' 


1 


Tr 


0 JjT 0 


V (t) dt = 


A 


Vnt, 


(j+l)To N 

T, COS (k„x - LO n t + </>„) = 
o JjTo n=1 


N 1 

y - [sin ( k n x - ui n (j + 1) T 0 + 4> n ) 

z ' UJ 


^ n—l 

sin (k n x - uj n jT 0 + 4> n )\ = 

N 


(3.2) = 2 


A 


y — sin f y n T 0 ) cos f k n x - 


VNTq 


n—l 


j + - ) T 0 + . 


This defines the Hamiltonian on one time interval; accordingly we can write the 
Hamiltonian of one interval as 


( 3 . 3 ) 


H ( o j) W = -l^2+Vo U) (*) 
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The global Hamiltonian corresponding to (2.1) of |T], which gives the zeroth order 
approximation to H (f), is generated by: 

(3.4) H 9 = H^\t) for jT 0 <t<(j + 1 )T 0 . 

Using this notation, a general time evolution can be written as the product of the 
interval propagators since Hjj is piecewise constant in time, 


( 3 . 5 ) u 0 (t) = e -^ ma *\t- jma *T 0 ) _ _ _ e -m^T 0 e -m^T 0 

where j m ax — 1 is the integer part of t/To The evolution in this order is 

(3.6) {x, t) = U 0 (i) ip (x, 0). 

The propagator satisfies 

(3.7) i^-U 0 (t)=(3H { 0 9) U 0 (t) 

corresponding to (2.10) of pQ. This propagator is numerically implemented and 
used to solve the time dependent equation of motion. The error in this order is 
bounded by /3? up to times of order T 0 , the reasoning behind this is shown in pQ, 
Eq. (2.49) there. The error in diagonalization of Hq' > is negligible (see discussion 
in Sec. 5) 


3.2. First order. The first order averaging is based upon a “peel-off” transfor¬ 
mation of the zero order. In such a transformation the next order Hamiltonian is 
constructed from the zero order in the following way: Let Hi (t) be defined as 

(3.8) Hr ( t ) = Uo 1 ( t ) [H (t) - H 0 9 (t)] U 0 ( t ) 


Hence Hi ( t ) is the Heisenberg dynamics of the full problem with Uq (t) dynamics 
peeled off. An important point to note is that the Laplacian term drops in Hi\ 
Hence Hi (t) is a bounded operator, for which the results of [1] directly apply. Its 
average in the j-th interval is 

_ ,1 rU+ 1 ) 

(3.9) H i =t Hl ^ dt 

1 0 JjT 0 

The Hi (t) dynamics is given by the propagator 

(3.10) Ui (:t ) = _ e -i0H[^T o 


We turn now to calculate H [^, by (2.49) of [T] it is of order Before diagonalizing 
this operator in order to calculate the time evolution, there are several steps needed 
to be taken. First we write explicitly Hi using integration by parts 

13-H) 

Hi U) = \r [ U+ } ° Uo\t) [H ( t ) - H% (t)] U 0 (t)dt = Hi W) + Hi U ’ U) + Hi (j ’ U1 \ 

1 0 JjTo 

where 


Hi 


(J, i) 




U 0 (t) 


(3.12) 


u + m 

t = jTo 
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(3.13) 

_ / ■ tt'i 1 /’(”+ 1 ) T 0 / F ) \ \ ft ' 

Sl = ~T 0 JnT dt \df U ° 1{t,) ) L l H ( s )- S o(s)]ds Uo(t'), 

and 

(3-14) 

, ttt'i I b ft ' _ / p ) \ 

Sl = ~nJ jT U ° {t,) i l H ( s )~ S o(s)]ds (t')dt'. 

We note that for any integer t/T 0 , 

rU+l) T o 

(3.15) / [H ( t') - H° {t')\ dt ' = 0 

JjT 0 

by construction. Therefore for such t = q and ( t ) can be simpli¬ 

fied. Moreover the expression in Hi^’ lir> is just the hermitian conjugate of Hi^’ ir> 
so we only need to analyze the second term: 

Hi ( WlQ)= jr J dt' J [Hj (s) - Ho (9) (a)] ds U 0 (t') 

i rti+i)T 0 / f) \ [ /’ f/ t \ 

= AT / dt'i wjUq 1 ^') / H 0 (s) ds - {£ - jTo) Ho { 3) u 0 (f) 

1 0 j jTo \<tt ) L J jTo 

To evaluate this operator we can use a recursive construction of states based on the 
zeroth order eigenstates of the averaged Hamiltonian. It is important to remember 
that this is actually a matrix. For convenience of notation we will define 

(3.17) Uo C t ) = . e -iPH^To e -iPH^T 0j 

For t that is an integer multiple of Tq. 

j 

(3.18) Uo (i) = J] Wj-j. 

f =o 

where 

(3.19) Wj(t) = e - i P B ° )r °. 

First we assume t/To is integer and then generalize for any arbitrary t. The eigen- 
values and eigenfunctions of Hq are 

(3.20) H& |^> = Ej |,■$) 

Here we use a base of size M that is assumed to be finite (in this work we take 
M = 64). Since /? is small the eigenfunctions are approximately eigenfunctions of 
a free particle, namely 

(3.21) < x\v k j >« e fkx . 

We will have M eigenvalues and eigenvectors, the largest value of k is 50. Between 
each pair of propagators W 3 . Mu_i we can insert the identity resolution in the 
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corresponding basis (3.201 
(3.22) 


M 




fe =1 


resulting in 

(3.23) U 0 (t) =W j ---W 1 W 0 = Wjl) ■ ■ ■ WihWolo 

As an example let us take just the first two terms 

M M 

WxhWoh = E W 1 1^ 1 ) (vM E W ° 1^0°) {<P\ 

fc 0 =i 


(3.24) 


k i=l 
fci | 


the brakets expressions (^i 1 1 Sfc 0 =i l^o”) are J us ^ scalars, and since ipj 
eigenfunctions of H 3 0 ( |3.20[ ) 

M M 

(3.25) WihWolo = E W ' l^ 1 ) (^i 1 I w o E l^o°) <¥>S°I 


k i = l 
M 


ko — 1 


M 


E e_l/3W 1^' 1 ) E l^o°) (a 

ko — l 


k i = l 


using the notation 
(3.26) 

this becomes 


a k 1 k 0 — {^l 1 \<p 0 °) 


M M 


(3.27) Wi/iWq/o = E E a kik 0 e-~ i ^ T ° El1 e~ i ^ T ° E °° l^ 1 ) <^g°| 

ko — 1 k i = l 


In the same way for the complete sequence of propagators (3.181 
M M 

Uo {t) = E ''' E afc j fe j-i Q:fe 3-i fe j-2 • • • 


ko — 1 kj = 1 


(3.28) 


and it satisfies (3.7). 

Here jT 0 < t < (j + 1) T 0 . 

The inverse takes the form 
(3.29) Uq 1 ^) = e ^ 0e ^)To. e *PH 0{ t- 3 T 0 ) 

M M 

= E-" E < fcl ^ lfc2 ^ (t - jTo)E °e^ (i - 1)£o ko fe °)(^ 


kj—1 k j — i — l 


By (3.7) or direct differentiation of ( |3.28[ ) 

r, M M M M 

( 3 - 30 ^ c/ °( < ) = -E ■ E E ■ E 

fen = l fc'=lfc0 = l fci = l 


fc 0 


e -m- jTo )E/ | *i } ^fca| 


are 
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and in a similar way one finds 
(3.31) 

p. MM 

s qr‘M• 

kj=l kj — i = l 

substitution of ( 3.30| ) and ( |3.28[ ) into ( |3.16[ ) leads to 


^koki^kik2 * ’ ’ ^kj — ikj^ 


•P T o E o # _ e i/9(*-jT 0 )S/ 


* H3EZ* |^°) (v7 


i?i W,II) (t) 


M MM M 

4 e-ee-e^ 


'j a k' 0 k' 1 ■ ■ ■ 0l k'-_ 1 k l j a kjkj-i ■ ■ ■ CKfcifco 


feo = l fc'=l feo = l %=1 


/'(.) +FT) k' k' fc' / , N A,' 

/ dt' e iPToE o e -^T° E o 0 . e i0(t-jTo)E/ e -ip{t'-jTo)E/ 

J 1T 0 

Wo°){vf 


[ H(s)ds-(f ~jT 0 )H 0 U) 

JjTo 

. MM M 

4 e-ee-e^ 


k?)wi= 


fco =1 fe )= lfe o=l fci=l 
"0'+l)T 0 


j a k' 0 k' 1 ■ ■ ■ OL k' j _ x k'. OL k j k j _ 1 ■ • ■ «fcife 0 


(3.32) 


'iT 0 

(b' 


dt' e -iP T o E o° ' < _ e iP{t-jT 0 ) E / e -i^(t-jT 0 )E j 

rt' 


H (s) ds-(t- jT 0 ) H 0 U) $ ) Wl°) <¥>o 


kn \ 


O'To 


o fe 0| 


Finally the full expression for the first order averaged Hamiltonian is, where (3.11) 
and (|3. 151 were used, 


(3.33) 

Hi (j) (i) 

(3.34) 


(t) = /r 1 w,n) (<) + ( 


T (1,11) 


H 


(*) 


®( 4 e-ee-e^ 

w - fc'-=ife 0 =. 

*CH-1)T 0 


M 


M M 


M 


*0 = ! fc'=lfc 0 = l fej=l 


Qfiu/ j,/... Olv.t b,/ 


&kjkj — i • • • 

r t' 


dt'e~ i0T ° E o° ... e if3{t ~ jTo)E i J e ~ il3< ' t - jTo)E j 


jT 0 


vf 


ljT 0 


dsH (s) - (t - jTo) F 0 0) 


fe ) 


fc n 


1^0°) <?5 


fco I 


The integral (3.34) is performed numerically as follows. The domain of integra¬ 
tion is derived into squares of size At x At. The integral is approximated by a sum 
of the values of the integrand of the middle points of the squares, multiplied by At 2 . 
The error in each term is of the order of At 2 . This can be improved substantially. 
Here it is not required since we can obtain the required precision in this simple way. 

At first sight this expression (3.34) might seem very complicated but can be 
understood quite easily; in fact what we have here is a matrix constructed from 
the sum of matrix products; the terms of the matrix involve only the products of 
eigenvalues and eigenvectors of H q . The benefit of calculating the propagator in 
this manner is that one only needs to diagonalize the Hamiltonian in each interval 
only once. As a result of (3.2) the averaged potential is of order yffi and a. k ^ k . are 
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elements of matrices that are almost diagonal as will be verified aposteriori in Sec. 
5. Consequently the error in the calculation of (3.34| is of the order 

(3.35) SI t = At 2 

By the general theory see (2.15) of [T| the first order propagator takes the form 

r ~— - (j) (j) 

(3.5) with Hq replaced by H\ ' namely with 

U 1 = e -iPH[ j) {t-jT 0 ) . e -i / 3J?( 1 'T 0e -i/3H< 0) T 0 

3.3. Normal form transformation. To improve the accuracy we implement the 
normal form transformation; we will use the form given in Eq. (3.10) of PQ. In our 
notation this becomes 

rt 


(3.36) 


t/i (t) = 1 + i0Ui (t) 


dt' 


H x (t!) - H[ 9) ( t') ) Ui (t). 


The manner in which we will simplify the above expression will be similar to the 
method used to calculate U\ and H\ given in Eq (3.10) and (20), splitting into 
intervals of length T 0 and using (3.151, replacing H 0 by Hi: 

(3.37) U\ (t) = l + ifiUi 1 


( j r 0'+ i) T ° 

E/ df 

Hi (t') - R[ j,) 

\ r=0 J J'T 0 

- 


IjTo L 


by definition of H\ 


U) 


Hi (; 1 0 + H[ o) (; 1 0 dt' ) Ui (t) . 
-(j'+l)T 0 


dt 1 


Ij'To 


Hi (t') - H\ 


O') 


= 0 


(3.38) Ui(t) = 1 + i/3Ui 


dt' 


tjTo 


Hi ( t') - H[ j) (t')] Ui ( t ) = 


i + mi 


1 jTo 


Hi(t')dt')-(t-jT 0 )H± 


U) 


Ui(t) 


and explicitly 
(3.39) 


U (t) = 1 + ipe^To ... e W'(t-i7b) 


U) ( 


1 jTo 


Hi(t')dt') -{t-jT 0 ) HI 


U) 


Ui(t) 


u ( t) = _ e -i(3H^T 0 ^ 

If also the normal form transformation is performed the error is of order /3 i . 

4. Iterative application and error analysis 

To reduce the error we introduce an iterative process. After the normal form 
transformation is performed we introduce a new Hamiltonian 


(4.1) 

where 

(4.2) 


1 




-1 


Hi-Hi 


(9) 


Ui 


Ui = U 0 UiU 
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The wave function at the new level satisfies a Schrodinger equation like (3.12) of jl] 
with Hi playing the role of A (t) there, and /3 is replaced by /3 3 / 2 . Now one starts 
from an equation like (2.13) with ft replaced by /3 3//2 and H by H\. At each step 
the effective value of ft is reduced (/? —► /3 3//2 ), and so is the error. The process is 
repeated until the bound on the error is satisfactory, as will be shown below. 

Assume the process repeated l times. The resulting approximation for the wave 
function is 


(4.3) (*) = Ui . Ui-iUrfiQ) 

where f/;/ is the propagator at the /' level of the hierarchy, corresponding to the 
Hamiltonian Hi. 

We turn now to estimate rigorously the errors using the multi-scale and averaging 
method, assuming /-levels of the hierarchy. 

With p small on a time interval of order ft? = Tq. Let us denote by t max the 
maximum time we want to simulate dynamics. After introducing a normal form 
transformation we eliminate the c/3$ error term on (3.10) of [I], and we get for the 
evolution 

(4.4) i/i(t) = U (t) ip (0) = U 0 (f) Ro ( t) ip (0) 
leading to 

(4.5) U(t) = U 0 (t)R 0 (t) for 0 < t < t m ax 

with Rq (t) = 1 + 0 for 0 < t < t max with U 0 (t) is given by (a product of) 

averaged dynamics followed by a normal form unitary transformation. Moreover 

(4.6) i-^=phh(t)R 0 (t) , 
with 


(4.7) 


|ffi(i)||=0(l). 


If we then use the same method of averaging on Rq (t) in Eq. (4.6) we get Ro(t) = 
Ui(t)Ri(t) and 


(4.8) 


U ( t ) = U 0 (t) Hr (t) Ri ( t ) 


where now U\{t) is the averaged approximate solution for Ro(t ), and 

(4.9) R 1 (t) = l + 0 (/?i) + O ((>3) 5 t ma ^j . 


Again after normal form transformation, the O 



correction drops and we have 


(4.10) U(t)=Uo(t)U 1 (t) + o([piy t max ^ . 

After /-such iterations, we get the exact solution 

(4.11) U (■ t) = U 0 (t) lh (*)••• Vi (t) + O (^)'i ram ) . 

So the convergence of the scheme is super exponentially fast, close to the Newton 
type iteration. 
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We denote by t max the maximum time we want to simulate our dynamics. The 
error in the MSA level l of the hierarchy denoted by e can be written as 

(4.12) e = p(-i) l t max . 

We can then invert this relation to obtain the number of desired hierarchy steps l 
for a given (3 and desired error e, 


(4.13) 


l > 


log ( 


log e —login 
log 0 




in terms of T max = /3t max (see Eq. 2.81 


(4.14) 


l > 


log( 


loge-logT n 

log (8 


+ 1 


log I 


The error in the integral (3.34) is given by ( |3.35[ ) and the error in the diagonalization 
of the averaged Hamiltonian is assumed to be small (see discussion in the end of 
Sec. 5). 


5. Numerical implementation 


The purpose of this section is to demonstrate that the multi-scale averaging 
method developed in [T] is superior to the standard split step method as it is much 
faster and the bound on the error can be estimated analytically. We will evolve the 
wave function for the model presented in Sec. [2] with the approximate evolution 
operator of the multiscale and averaging (MSA) method and compare the results 
to the ones found using the standard split step method. In the zeroth order we 
evolve the wave function with the help of ( |3.6| ), with Uq calculated by (3.28). 

The diagonalization (3.20) can also be performed once and can be done in parallel 
for the various time intervals. In particular the first order (3.10) requires to di- 
agonlize H\ 3> that in turn is given by the diagonalized Hq ^ given by (3.3). This 
enables to compute the ,k 0 of (3.26). To obtain the first order MSA approximate 
dynamics, we use the evolution operator as given by (4.3). 

The results are compared to the ones found with the help of the split step method. 
In this method the wave function is propagated keeping only the kinetic energy or 
the potential energy in small steps of size St. The value of the potential is taken 
in the center of the time interval. The choice of a time step St is crucial. The way 
to test the convergence of the scheme is by running the dynamics up to a point t 
using a time step St, running the dynamics up to the same point t only using a new 
time step St' = If the wave function is the same within some fixed accuracy 
then the scheme is assumed to be convergent. The accuracy is defined as 


(5.1) S a (t)= J dx\il>[ St) (x,t)-ip { 2 2 \x,t)\ 2 . 

Where T is the domain in space where the wave function is defined. If the error 
S a (t) is not small then one needs to continue adjusting until one converges. Listed 
in Table [T] are some values of the small parameter j3 and the time scale it dictates. 
For longer time scale a smaller and a more refined time step is needed in order 
to converge the split step method. For the values listed in Table [T] an accuracy 
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of S a = 10 7 was chosen as a convergence criterion, i.e if the two wave functions 
obtained at the same time t with different time steps St and differed by less 
then S a = 10 , the algorithm is considered converged and an appropriate choice 

of St is obtained. The * marks the fact that the required time was too long for the 
standard split-step computation to converge. 

To demonstrate the accuracy of the MSA method, we denote by i[) a (x,t) the wave 
function computed by the MSA method and by ijj s (x,t) the one found by the split 
step method and compute the deviation 


(5.2) A(t) = J dx\il> s (x,t) - if} 0 {x,t)^ 

The initial wave function in all our computations is 

x 2 

(5.3) i>(x,t = 0) = 


where J\f is the normalization constant. In Fig. |6.1| the comparison between the 
wave functions ip s ( x , t) and ijj a (x, t ) is obtained by evolution starting from the 


initial wave function (5.3). It is presented for an arbitrary realization of the random 


potential (2.11). The difference is very small and it will be calculated in what 


follows. The results are presented in Fig |6.2| where we plot 
(5-4) A (t) = (A (£)) av . 

The average is over 40 realizations of the potential V (x,t) of (11). The averaging 
is performed as follows: A (t) is calculated for a specific realization and the average 
is taken so that the (f>i are distributed uniformly in the interval [—7r, 7r], while Vn ^ 
and k n are distributed uniformly in the intervals [v mm = — 15,v max = 15] and 
[k mm = — 20, k max = 20] respectively (see Sec. 2 paragraph following Eq. 5 ). ui n 
is calculated by (|2.5|). We take /3 = 1 ■ 10~ 3 and /3 = 1 • 10 -4 while <r° = 1 resulting 


in the initial momentum standard deviation ct 3 = 0.5. In our basis Eqs. (3.20) and 


(3.21) for \ifij >, the largest value of the momentum is 


= 80, therefore 


(5.5) 


k 




The largest value of |a;| is 10 therefore cr° is much smaller then the largest value of 
|x|, namely 10. Note that A (t) is much smaller then S a and e, indicating that the 
results are much more accurate then expected from the theoretical bounds. This 
may be specific to the potential we used. A similar situation was observed in the 
appendix of jl]. 

To demonstrate the efficiency of the calculation we compare the computer time 
Tcomp required to perform the numerical time propagation Fig. [2] up to times 
t-max = 8000, we compare the results for the multi-scale and averaging split steps 
methods with the same precision S a = e = 10~ 7 we choose a time t max = T 0 ■ i 
such that i is the smallest integer satisfying t max > t max . the lowest hierarchy l m i n 
required for the calculation l is used. The results are summarized in table [2] and 
plotted in Fig |6.3| The calculations were performed on two computational nodes 
each composed of a 2.4 Ghz Intel Xeon processors. 
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p 

V 

rp 2 

1 0 

1 0 

St for Tq 

St for Tq 

1•10" 2 

100 

100 

10000 

0.01 

0.0009 

1•10" 3 

1000 

1000 

1000000 

0.01 

0.0001 

5•10- 4 

2000 

2000 

400000 

0.009 

* 

1•10 -4 

10000 

10000 

100000000 

0.009 

* 


Table 1. split step time steps for different values of the small 
parameter. The accuracy is S a = 10~ 7 


It turns out that for the problem we studied, the results we obtained are probably 
much better then the error estimate (4.12). To see this we present in table [3] the 
difference 


(5.6) 


A u+1 (t) = J dxl^x.U 


c) tmax) \ 


indicating the order of magnitude of the error as well as the bound (4.12) for l = 3 
and l = 4 for different values of /3 and t max = Tq = 4 Indeed this is the case. 


The error in the calculation of the integrand (3.34) is given by (3.35). The reason 
is that the matrix consisting of the ot ki , kj is almost a unit matrix. For all of our 
calculations we verified that 


(5.7) 

M—64 

\o-ki,kj < 10 -23 

while 

kj y^ki 


^ \a kukj -1| < 1(T 20 . 

kj=ki 

The error in the cliagonalization of the averaged Hamiltonian Hj 9 ' 1 is of the order 
of 10 -50 . The diagonalization is performed by the lancos algorithm | B41 . 


P 

5•10" 5 

1•10 -4 

3•10~ 3 

1•10" 2 

3-10" 2 

1•10" 1 

T 0 

141.42 

100 

18.25 

10 

5.77 

3.162 

i 

57 

81 

439 

801 

1386 

2530 

^max 

8061.02 

8100 

8015 

8010 

8002.07 

8000.56 

Imin 

3 

3 

4 

5 

5 

6 

rpav 
^ comp 

93 

63 

50 

25 

8 

4 

fpss 
^ comp 

1140 

1080 

780 

580 

420 

380 

comp 

T ss 

comp 

0.081 

0.058 

0.064 

0.043 

0.019 

0.011 


Table 2. The comparison between the computational time T°^ rnp 
(in minutes) using the multi-scale averaging method and the com¬ 
putational time T** mp using the split step method. Averaging over 
40 realizations similar to the ones used in Fig. 2 was performed. 

































NUMERICAL IMPLEMENTATION OF THE MULTISCALE AND AVERAGING METHODS FOR QUASI PERIODIC SYSTEMS4 


p 

5 • 10~ 5 

1 ■ 10 -4 

3 • 10 -a 

1 • 10~ 2 

3 • 10 -2 

1 • 

io - 1 

CO 

<1 

6 • 10 -11 

1 . 04 - 10" 10 

1 . 44 - 10" 10 

1.81 • 10" 10 

2 . 17 - 10" 12 

2 . 5 - 

10-ro 

Q=3 

75 • 10 ~* 

3.88 • nr 7 

2.02 • 10~ e 

1.78 • 10 b 

2 . 42 - 10" 4 

2.21 

■ 10" a 

e;=4 

6 . 5 - 10" 10 

1.08 • 10~ 9 

1.81 ■ 10“ 9 

7.5 ■ 10 -8 

6.5 ■ 10 -7 

8.66 

■ 10" 5 


Table 3. The difference (5.6) and the bound (4.12) for various 
values of /3 and t max = T 0 2 = g. 


6. Spreading in k-space 


Multiscale and averaging (MSA) enables us to calculate very accurately the 
spreading in momentum space over a very long time. For this purpose we evolve 


the wave function ip ( x,t ) starting from (5.3) with a x = 1 for the potential (2.11). 


The wave function is used to calculate the variance of the momentum 


Var k{t) = j ip* (k,t) (k — k) 2 ip (k,t) dk 


k= ip* (k,t) kip (k,t) dk 


( 6 . 1 ) 

where 

( 6 . 2 ) 


and ip ( k , t) is the Fourier transform of ip ( x , t). Then we calculate the spread of the 
momentum relative to initial one 

<“> AVa 4^ V "^r (0) 

This calculation is performed for each realization of the random potential. Then 


average over 40 realizations of the random potential was performed as in (5.4), 
namely we calculate 


(6.4) 


(t) = (AVai{ n) (t) 


It is plotted as a function of r = pit in Fig. |6.4[ Because of the smallness of A„ 
we conclude that the replacement of the Hamiltonian by a finite matrix does not 


effect the result (see Eq. (5.5)). The plot is smoothed by a verag ing over intervals 
of length At = 10 2 , leading to the results presented in Fig. 6A the calculation is 
repeated for several values of /3, and in Fig. |6.6|the results are fitted to the formula 


(6.5) A„(r) = C(/3)r a ^. 
From the plot it is reasonable to extrapolate 

( 6 . 6 ) 


lim C = lim a = 0. 

/3—>0 /3->0 


We conclude that in the limit /3 —» 0 the spreading stops. This is the limit where 

(r) i— 

the velocity is much larger then the Chirikov resonant velocities of Eq. (2.5). 


(r) - 

It leads us to the conjecture that if in one dimension the v n ; are bounded, the 
kinetic energy cannot grow to infinity, in-spite of the driving. 

In Fig. 6.7 the classical and quantum results are compared. The quantum 


results were computed as the ones for Fig. |6.4| while for the classical results 60 initial 
conditions were also chosen at random from a Gaussian distribution corresponding 


to the initial quantum wave function ip(x,0) of (5.3). Both classical and quantum 





























NUMERICAL IMPLEMENTATION OF THE MULTISCALE AND AVERAGING METHODS FOR QUASI PERIODIC SYSTEMS5 


results were averaged over 40 realizations of the random potential. We note that 
both classical and quantum spreading in momentum stops. The classical spreading 
stops at an earlier stage. 
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Figure 6.1. The wave function at time t = Tn = -U obtained 

VP 

from the averaging method with precision of e = 10“' and the 
split step method with accuracy of S a = 10“' for /3 = 10“ 3 and 
cr° = 1.0. Two cases of multi-scale and averaging are presented (a) 


The fourth level of the hierarchy (l = 4, Eq. 4.12) (b) The fifth 


level of the hierarchy (l = 5, Eq. 4.12). 
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Figure 6.2. The average difference A (t) (Eqs. 5.4 (5.2)) between 


the split step and averaging method as calculated for two different 
values of (a) (3 = 10 -3 and (b) p = 1CT 4 while cr° = 1.0. Each 
subfigure presents both calculations for the third and fourth level 
in the hierarchy (l = 3,4 Eq. 


e = 8 a = 10 


-7 


4.12). The precision required is 
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Figure 6.3. Computational benchmark (run time) as a function 
of /3 as presented in Table [ 2 ] (a) Run times, obtained for split step 
method (T^ mp ) and for the multiple scale and averaging ( T'comp) 
method. The precision required was e = 6 a = 10 -7 and the values 


of l are presented in table 

12 


(b) The ratio 


presented in table 
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Figure 6.4. The average spread of the kinetic energy A„ of Eq. 
(6.41 of a wave function obtained for the multi-scale and averaging 
method with (a) /3 = 5 • 10~ 3 and l = 9 (b) (3 = 1 • 10~ 4 and l = 8 
. For both cases e = 10 -10 . Note that the hierarchy used here is 
much higher then required for the precision e. 
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Figure 6.5. A„ (r) smoothened over time intervals At = 10 2 
as a function of r on (a) regular scale (b) logarithmic scale, for 
/3 = 10 -4 , e = 10~ 10 and 1 = 6. 
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Figure 6.6. The parameters of the fit (6.5) as a function of /3 
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Figure 6.7. The average spreading of the quantum and classical 
kinetic energy. The quantum results were found with the multi¬ 
scale and averaging method and the classical result is obtained via 
standard R.unge Kutta integration with integratio n th reshold of 
10 -10 . The parameters used are the same as in Fig. 6.4 but l = 6 
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